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ABSTRACT 

Quadratic methods with heuristic weighting (e.g. pseudo-C; or correlation function 
methods) represent an efficient way to estimate power spectra of the cosmic microwave 
, background (CMB) anisotropies and their polarization. We construct the sample co- 

' variance properties of such estimators for CMB polarization, and develop semi-analytic 

techniques to approximate the pseudo-C/ sample covariance matrices at high Legendre 
multipoles, taking account of the geometric effects of mode coupling and the mixing 
between the electric (E) and magnetic (B) polarization that arise for observations 
covering only part of the sky. The E-B mixing ultimately limits the applicability of 
heuristically-weighted quadratic methods to searches for the gravitational-wave signal 
in the large-angle -B-mode polarization, even for methods that can recover (exactly) 
^ (— | unbiased estimates of the _B-mode power. We show that for surveys covering one or 

two per cent of the sky, the contribution of _E-mode power to the covariance of the 
q , recovered i?-mode power spectrum typically limits the tensor-to-scalar ratio that can 

' be probed with such methods to ~ 0.05. 

^ | Key words: cosmic microwave background - methods: analytical: - methods: nu- 

merical. 



1 INTRODUCTION 

With the recent results from the Degree Angular Scale Interferometer (DASI; Kovac et al. 2002; Leitch et al. 2004), the Cosmic 
Anisotropy Mapper (CAPMAP; Barkats et al. 2005) and the Cosmic Background Imager (CBI; Readhead et al. 2004), and the 
imminent two-year results from the Wilkinson Microwave Anisotropy Probe (WMAP; Kogut et al. 2003), the characterization 
of the polarization of the cosmic microwave background (CMB) polarization is gathering pace. The ultimate goal of such 
experiments is to constrain parameters describing the cosmological model and its fluctuations, and hence test models for the 
generation of these fluctuations, such as inflation. In Gaussian theories, all of the cosmological information encoded in the 
CMB anisotropies and their polarization is contained in their power spectra. For this reason, estimating CMB power spectra 
is a useful and efficient form of data compression in analysis pipelines (Bond, Jaffe & Knox 2000). However, estimating the 
spectra is not enough: it is also essential to have an accurate covariance matrix, and preferably an accurate model for the 
shape of the likelihood function (see e.g. Bond et al. 2000; Efstathiou 2004). Modelling the covariance matrix for a particularly 
efficient class of polarization power spectra estimators - heuristically-weighted quadratic estimators, or pseudo-C;s - is the 
subject of this paper. 

The polarization field can be decomposed uniquely on the full sky into a gradient part, with electric (E) parity, and a 
curl part with magnetic (B) parity (Kamionkowski, Kosowsky & Stebbins 1997; Seljak & Zaldarriaga 1997). The cosmological 
importance of this decomposition is that linear density perturbations cannot produce magnetic polarization, and so, in standard 
scenarios and at linear order, large-angle B-mode polarization is produced solely by gravitational waves. Second-order effects, 
most notably weak gravitational lensing by large-scale structure, can generate B modes from primary E modes (Zaldarriaga 
& Seljak 1998). The lens-induced B modes have a blue spectrum, and should dominate any primary signal from gravitational 
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waves on sub-degree scales. On larger scales, they may still dominate if the gravitational wave amplitude is sufficiently low. 
Either way, the variance of the B modes is expected to be much smaller than that of the E modes. Given the different parity 
properties of the E and B fields, CMB polarization in a statistically-isotropic and parity-respecting universe is characterized 
by two power spectra, Cf and C ; s , and also the cross-power Cf E between E and the temperature anisotropies. 

A number of different methods have been developed to estimate these CMB polarization power spectra. The opti- 
mal method, brute-force maximum likelihood (e.g. Gorski f994; Bond, Jaffe & Knox f998), and its quadratic approxima- 
tion (Tegmark & de Oliveira-Costa 2002), should produce the most accurate (minimum- variance) estimates. However, these 
methods require the inversion of iV p i x x 7V p i x matrices, where iV p i x is the number of pixels in the map, which is prohibitively 
slow for the analysis of mega-pixel maps. Given these limitations, a number of fast alternative methods have been sug- 
gested (Hansen & Gorski 2003; Kogut et al. 2003; Chon et al. 2004) building on earlier work on the analysis of the temperature 
anisotropies (Szapudi, Prunet & Colombi 200f; Wandelt, Hivon & Gorski 200f ; Hivon et al. 2002; Hansen, Gorski & Hivon 

2002) . These fast methods are closely related, and all involve compressing heuristically- weighted maps of the Stokes param- 
eters into a set of pseudo-C;s making use of spherical transforms - a step that requires only O(N^) operations with fast 
Fourier transform methods on iso-latitude pixelisations such as the widely-used HEALPix (Gorski et al. 2004). The methods 
differ in the details of how they transform the pseudo-C;s to power spectrum estimates. In the absence of instrument noise, 
the mean of the pseudo-C;s is linear in the true power spectra, so, for observations over a large enough part of the sky, 
unbiased power spectrum estimates can be obtained by applying the inverse operation to the pseudo-C;s (Kogut et al. 2003). 
If the observed part of the sky is sufficiently small, direct inversion is not possible and some form of regularisation is required. 
Hansen & Gorski (2003) advocate using likelihood methods, while Chon et al. (2004) make use of correlation functions as 
an intermediate data product. A notable feature of the latter approach is that the inversion is constructed so that electric 
and magnetic polarization power is not mixed in the mean. Pseudo-C; methods have now been applied to a number of real 
data sets, both for temperature (Netterfield et al. 2002; Hinshaw et al. 2003) and polarization analyses (Kogut et al. 2003; 
Ponthieu et al. 2005). 

Despite the attention that pseudo-C; techniques have received, relatively little work has been done on trying to understand 
their error properties. This paper aims to fill that gap by providing the basis for a reliable, semi-analytic error analysis technique 
for CMB polarization power spectrum estimates obtained with pseudo-Ci methods. This extends earlier analytic work that 
discusses temperature anisotropies only (Hinshaw et al. 2003; Efstathiou 2004), and the recent thesis work by one of us (Chon 

2003) where both temperature and polarization are considered. A significant complication that arises for polarized data is 
the mixing of E and B modes for observations covering only part of the sky (Lewis, Challinor & Turok 2002; Bunn et al. 
2003). We develop approximations to the pseudo-C"; sample covariances, taking this mixing into account at leading order for 
maps with smoothly-apodized edges. Our approximations are accurate in the regime where the mixing is perturbative, i.e. for 
I >• 'max, where Z max is the effective band-limit of the function used to weight the Stokes parameters. It is difficult to construct 
good approximations for the covariance on large scales, I <C i m a X - This is partly because of the increased importance of E-B 
mixing there, but also because the polarization power spectra do not go smoothly to zero on large scales in the presence of 
reionization. (The behaviour of the temperature-anisotropy power spectrum on large scales produces a similar effect for the 
temperature pseudo-Ci covariance.) Following Efstathiou (2004), the \ow-l sector of the covariance matrices can be constructed 
directly with exact methods, and used to overwrite the inaccurate approximation. Being able to evaluate accurate covariances 
efficiently is important in any application where multiple evaluations are required. Examples include parameter estimation, 
where the covariance generally depends on position in parameter space, and optimising the design and analysis of surveys, 
where one might wish to assess a large number of weighting schemes. 

Throughout, we ignore the effects of instrument noise so we deal only with the sample covariance of the C\. This is to 
highlight the geometric effects of heuristic weighting that are peculiar to polarization analysis, particularly E-B mixing. An 
obvious extension of this work that will be required for application to real data is to include the contribution from instrument 
noise. For uncorrelated noise this should be straightforward, but in a more realistic set-up, one will probably have to resort 
to Monte-Carlo simulations for the noise contribution. 

As noted above, one of the future goals of CMB polarization observations is to detect (or constrain) the stochastic 
gravitational wave background via its signature in large-angle B-mode polarization (Seljak & Zaldarriaga 1997; Kamionkowski 
et al. 1997). Magnetic polarization allows much smaller amplitude gravitational wave backgrounds to be detected than with 
temperature or electric polarization alone, since, in the latter case, the gravitational-wave signal is swamped by the cosmic 
variance of the signal from the dominant density perturbations. Methods to separate out exactly pure magnetic modes from 
observations with incomplete sky coverage exist at the map level (Lewis et al. 2002; Bunn et al. 2003; Lewis 2003). Such 
methods have the desirable property that quadratic power spectrum estimates formed from the pure-B modes have no cosmic 
variance if the B-mode power is zero. Were it not for gravitational lensing, separating B modes at the level of the map would 
thus allow the detection of arbitrarily low amplitude gravitational wave backgrounds in the absence of instrument noise. In 
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practice, imperfections in removing the lensing signal limit the amplitude to tensor-to-scalar ratios 1 r~3x 10 -5 (Kesden, 
Cooray & Kamionkowski 2002; Knox & Song 2002; see Seljak & Hirata 2004 for refinements from more optimal lensing- 
reconstruction methods). Heuristically- weighted pseudo-C; methods do not share the property that the cosmic variance of 
the estimated Cf is immune to _E-mode power, even if we take care to construct unbiased estimates so that e.g. the mean 
of our Cf estimator decouples from Cf (Chon et al. 2004). We investigate this issue further here, exploring the limits of 
detectability for r that arise solely from cosmic variance of the _E-mode polarization leaking into the variance of quadratic 
estimates of Cf . 

This paper is organised as follows. Section 2 reviews the properties of the polarization pseudo-C;s and establishes our 
notation. Following a brief discussion of the exact pseudo-C; covariance matrices in Section 3, we develop analytic approxi- 
mations to them that can be rapidly computed in Section 4. We begin with a toy model, based on a Gaussian weight function 
applied to a small survey area in the flat-sky limit. (Some of the more involved manipulations are relegated to Appendix A.) 
This instructive example illustrates the key issues that are specific to the analysis of polarization, i.e. the geometric mixing 
of E and B modes. We then move on to develop approximations for arbitrary (smooth) weight functions, and test these 
against the exact covariances in cases where we can easily compute the latter (in particular, for azimuthally-symmetric weight 
functions). Simple rules-of-thumb for band-power variances are also derived in Section 4.4 as limiting cases of the more general 
expressions; some of the manipulations of the 3j symbols required for these derivations are relegated to Appendix B. Finally, 
in Section 5 we consider the implications of our results for the detection of gravitational waves with B-mode polarization 
estimated with pseudo-C; methods. 

Throughout, we illustrate our results with the best-fitting power- law concordance ACDM model to the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) and 2dF Galaxy Redshift Survey data (Spergel et al. 2003; Table 7), but with the 
tensor-to-scalar ratio r = 0.15 corresponding to <f> 2 inflation. We include the effects of gravitational lensing in the B-mode 
power spectrum, but ignore effects due to its non-Gaussianity. Theoretical power spectra, pseudo-C;s and estimated power 
spectra in all figures include the smoothing effect of a Gaussian beam of ten-arcmin FWHM. 



2 PROPERTIES OF POLARIZATION PSEUDO-ClS 

We define Stokes parameters along the line of sight n using the basis vectors of a spherical polar coordinate system. We take 
the local x axis to be and the y axis — <p; these form a right-handed triad with the radiation propagation direction — n which 
defines the local z axis. The complex polarization P = Q + ill is then spin -2 (Newman & Penrose 1966) and can properly be 
expanded in terms of spin-2 spherical harmonics (see e.g. Lewis et al. 2002 for a review): 

(Q±iU)(n) =Y.(E^^iB lm ) T2 Y lm {n). (1) 

lm 

Reality of the Stokes parameters ensures that _E ; * m = (— l) m _E;_ m and similarly for B; m , while under parity transformations 
E lm — > (— 1)' Ei m (electric parity) but B lm — > (— l) l+1 B lm (magnetic parity). The polarization power spectra are defined in 
statistically-isotropic and parity-respecting models by 

(Ei m Eii m i) = C; 8iii5 mrn i, (Bi m Bii m i) = C l 5 H iS mm i. (2) 

The E-B power vanishes (and B is uncorrelated with the temperature anisotropics) by parity. 

In pseudo-C; methods, one adopts a heuristic weighting of the complex polarization with a scalar function w(n). The 
weight is necessarily zero in those regions that are not observed, but typically it is desirable to excise a larger region to remove 
further regions of high foreground emission. One can attempt to choose w(n) to improve the accuracy of subsequent power 
spectrum estimates. For example, on small scales, where the instrument noise dominates the signal, the optimal weighting 
is essentially inverse-noise- variance weighting (Efstathiou 2004), while on large scales uniform weighting is preferable. In this 
paper we shall assume that the weight is chosen to be real and to go smoothly to zero as excised regions are approached. 
The former restriction ensures that the weighting does not alter the polarization direction, while the latter makes w(n) 
approximately band-limited to some Z m ax- Such edge apodization allows us to deal with the effects of the resulting E-B 
mixing in the pseudo-C; covariance matrices in a perturbative manner; see Section 4. The problem of E-B mixing without 
edge apodization, i.e. taking w(n) everywhere one or zero, is discussed extensively by Lewis et al. (2002) and Lewis (2003) in 
the context of isolating B modes at the map level. 

Following Hansen & Gorski (2003), we define pseudo-multipoles Ei m and B lm to be the E and B-mode multipoles of the 
weighted polarization map, so that 



1 Our definition of the tensor-to-scalar ratio follows Martin & Schwarz (2000), i.e. it is the ratio of the primordial gravitational wave 
and curvature power spectra. 
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E lm ±iB lm = £ ±2l(im)(im)>(E(i m y ±iB(i m y), (3) 

(Im)' 

where ±2l(im)(lm)' are the Hermitian coupling matrices 

±2/(im)(im)' = dnw(n) ± 2Y (lm y(n) ±2 Y* m (n). (4) 



Using s Yj^ = {-l) s+m -sYi- m , and the reality of w(n), we have ±2l^ m )(imy = (-l) m+m T2 /(j-m)(i-m)' which ensures that 
E* m = (— l) m Ei- m and a similar result for Bi m . If we expand the weight function in terms of its spherical multipoles, 
w{n) = J2im w imYi m (n) , the coupling matrices can be evaluated in terms of Wigner-3j symbols: 



±2'(lm)(Im)' 



LM 



Vf n ™ / (2L + l)(2i + l)(2Z' + l) / / /' L W I /' L , 

WLM V 4^ m -m' -M ( T 2 ±2 J" (5) 



From equation (3) we can separate the E and B pseudo-multipoles as 

Elm = £ (+Ilm(lm.yE(lm.y +i—I(lm)(lmyB(l m y) , (6) 
(im)' 

-Bim = £ [+I{lm)(lm)'B(l m y — i-I(lm)(lm)' E(lm)') i (7) 
(im)' 

where we define the Hermitian ±I(i m )(im)' to be 

+-f(im)(im)' = 2 (+2^(im)(Im)' + -2-fym)(!m)' ) , -^(Im)(Im)' = ^ (+ 2 ^(!m)(!m)' _ -2^im(im)')- (8) 

Given the -E (m and -B (m , we form their pseudo-CiS according to 



Ci =27 



^^l^iml 2 , Cf ^^El^l 2 ' ( 9 ) 



We could also form Cf B = £ m Ei m Bf m /(2l + 1), but, as we show shortly, this vanishes in the mean in the absence of parity 
violations. However, it may provide a useful test on parity-violating physics or unaccounted-for systematics, and so future 
data should, of course, be validated with Cf B . 

In the absence of instrument noise, the mean pseudo-Cis (i.e. averaged over CMB realisations) are related to the true 
power spectra by convolutions. We find (Hansen & Gorski 2003) 

<Cf ) = E( P »' C ^ + M ^C?), <<5f > = £(M H ,C,? + PwCF), (10) 
c c 

where 



Pw 



21 



21 



^ £ M(im)(im)'| 2 = £( 2L + l)w L [l + (-1)*] ( 2 o) , (ID 

mm' L \ / 

^ T E|-/(im, ( imv| 2 = ^E(2L+lV i [l-(-l) K ]( j 2 2 o) 2 ' ^ 

mm' L \ / 



Here if = Z + Z' + L and we have introduced the power spectrum of the weight function wi = Yl m \ w im\ 2 / {21 + 1). A non-zero 
matrix M u i biases the pseudo-CiS by mixing e.g. i5-mode power into C B . For observations with uniform weight over the full 
sky, w(n) = constant and M H i = since then only wo is non-zero and the 3j symbol forces I = I' and hence K to be even. The 
sum of P H i and M H i is approximately equal to the matrix that relates the temperature CiS to the mean pseudo-CiS (Hivon et 
al. 2002) at high I. Their sum differs from the temperature case only by the presence of non-zero azimuthal quantum numbers 
±2 in the 3j symbol. The mean of Cf B evaluates to 

<Cf B >=£( P »'- M »')^ S , (13) 
r 

and so vanishes if parity is preserved in the mean {Cf B = 0). 

An example of the mean pseudo-CiS for observations over a circular region of radius 15° is shown in Fig. 1. The weight 
function is uniform inside a circle of radius 10°, but the remaining 5° is apodized as cos 2 [36(0 — 10°)] so that w(n) falls smoothly 
to zero at the edge of the observed region. For this weight function the matrix Phi is roughly two orders of magnitude larger 
than Mui, and both are narrow compared to features in Cf and C B except for the reionization signal on the largest scales. 
On small scales the mean pseudo-CiS can be approximated by removing Cf and Cf from the convolutions in equation (10) 
to give 
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Figure 1. Left: power spectra for B (top) and E (bottom; solid lines) compared to the mean pscudo-C;s (dotted lines) and the recovered 
power spectra (dashed lines). Right: the bottom panel shows representative window functions P ;i / (solid lines) and M ;i / (dotted lines) 
that give the mean pscudo-C;s on convolving with the true C;s; the top panel shows the pseudo-inverses P^, 1 (solid lines) and M~} 
(dotted lines) that when convolved with the mean pseudo-C;S remove the effect of E-B mixing. Note that M ;; / and M~} have been 
multiplied by a factor of 100 for clarity. The weight function applied to the map is uniform inside a circle of 10° radius, with cosine 
apodization out to 15°. To obtain the pseudo-inverses, a Gaussian apodization of 4° HWHM is applied to the correlation functions. 



<cf > « cf J2 p iv + cf J2 M »' . <cf > w c * E p »' + c ? E M »' • ( 14 ) 
i' i' i' i' 

For Cf the mixing term described by M u i is essentially negligible, and the mean of Cf is approximately a scaled version of 
the underlying power spectrum. The relevant scaling is discussed below. For Cf mixing is important, despite Pw dominating 
over Mui, since Cf 3> Cf ' . This additional contribution to (Cf) from geometric mixing of E and B modes can be clearly 
seen in Fig. 1, and, as expected, looks like a scaled version of the S-mode power spectrum. 
The normalisations J2i> Pw an d J2f evaluate to 

Y^Pw+M u , = £^±I Wi=w ( 2 >/ sky , (15) 

l' L 77 

„ t-^2L + 1 ( L(L + 1) (L + 2)\ (/-2)!\ , . 



4tt V 1(1 + 1) (L- 2)1 (1 + 2) 

where w^ 2 -*/sk y is defined by Airw^ / a k y = f w l (n) Ah, and / s k y is the fraction of sky with non-zero weighting. Equation (15) 
follows from the expansions of Put an d M H i in equations (11) and (12), and the orthogonality of the 3j symbols. It is the 
same normalisation as for the temperature anisotropies. Equation (16) can be established following the method used to derive 
equation (66) in Chon et al. (2004); see also Appendix B. Note that for / much greater than the assumed band-limit of w(h), 
the normalisation J2i> is suppressed by a factor 0(i max /i ) compared to Pw ■ This is consistent with the expectation 
that E-B mixing is suppressed at high I (Lewis et al. 2002; Bunn et al. 2003). We can gain further insight into this result by 
noting that [for smooth w(ii)] 

^(2L+l)L(L + l)w L = /"|Vw M | 2 dn, (17) 

E( 2L + 1 ) }r Z % WL = I «vV 2 (V 2 + 2) WM dn, (18) 

where w Bi i(n) = m wlmYlm^) is w(ii) smoothed with a top-hat (/-space) filter. For band-limited weight functions, 

at I 3> i max the smoothing has little effect and so to leading order 
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^ M »-lJwTT) dA {l>>l - y (19) 

At a given scale I, the relative importance of mixing of E- and B-mode power in the pseudo-Cjs is thus controlled by l/l 2 
times the ratio of the average squared gradient of the weight function to the average square weight function. 



2.1 Inversion of the pseudo-C;s 

For the purpose of parameter estimation, it is not strictly necessary to attempt to invert the pseudo-CjS to (unbiased) power 
spectrum estimates. Forming the pseudo-C;s can be regarded as a mildly-lossy form of data compression of the maps of 
Stokes parameters, and inferences about cosmological models can be made directly from them if their covariance - or, better, 
their full sampling distribution - can be calculated (Wandelt et al. 2001). The calculation of the signal contribution to these 
covariances is the main topic of this paper. However, for presentation purposes it is useful to be able to plot quantities that are 
simply related to the underlying power spectra. At the least this requires some approximate deconvolution of the geometric 
effects of the weight function from the pseudo-Czs. 

For observations covering enough of the sky 2 , the matrices P H i and M H i will be invertible. In this case we easily obtain 
unbiased estimates of Cf and Cf from 

Cf + Cf = £(P + M)-}(Cf + Cf), Cf - Cf = J2(P ~ M)w(Cf - Cf). (20) 
c c 

Alternatively, we can form unbiased estimates of the correlation functions of the Stokes parameters from the polarization 
pseudo-C;s (or directly from the maps), and then invert these with the appropriate Legendre transforms (Chon et al. 2004). 
For observations over smaller parts of the sky, P H i and M H / will not be invertible. In this case we do not have the spectral 
resolution to obtain unbiased estimates of the C;s at every I. Various approaches have been suggested to regularise the inversion 
including the use of band-powers (Hivon et al. 2002; Hansen et al. 2002), or apodized Legendre transforms of the (incomplete) 
correlation functions (Szapudi et al. 2001; Chon et al. 2004). The latter method provides estimates of the Cis convolved with 
a window function that depends only on the choice of apodizing function that is applied to the correlation function. Chon et 
al. (2004), building on earlier work by Crittenden et al. (2002), showed how to generalise the correlation function method to 
ensure that the estimated polarization power spectra did not mix E and B power in the mean. This is particularly useful for 
visualisation of the B-mode power spectrum which can easily be swamped by the dominant E modes. We do not reproduce 
the method of Chon et al. (2004) here, but note only that it provides a convenient way of constructing pseudo-inverses P^, 1 
and M7,, to Pw and M ;i ' respectively with the properties that 

Y,{PwPin» + M-}M VV ,) = - 2 K U „, Y,( p w M vi" + M^Pi>v) = 0- (21) 

V V 
These relations ensure that our estimates of the power spectra, 

Cf = E ( P w G v + Ma'Gf), Cf = J2 i P u>Cf + Mu'cf), (22) 
v v 

satisfy (Cf) = -iK u iCf and the equivalent relation for Cf . The normalised window function -iKiv = -2-fQj' / J2l -iKil, 
where -iKiy is given by (Chon et al. 2004) 

11' A- 1 r i 

- 2 K W = — — / /(/?)<£ -2(M-2G8)d cos /3. (23) 



Here, d l mn ([3) are the reduced Wigner rotation matrices and /(/?) is the apodizing function that is applied to the correlation 
functions. An example of the pseudo-inverses obtained via this route is shown in Fig. 1, along with the mean of the recovered 
C;s. The correlation functions have been apodized with a Gaussian with HWHM of 4°. The apodizing function was deliberately 
chosen to be narrow to ensure that the resulting pseudo-inverses P^, 1 and M~} are well localised. Given that the pseudo- 
inverses are significantly broader than P u i and M u i , the window function -2K H i (not shown in the figure) inherits the width 
of the pseudo-inverses. For this reason the spectral resolution of the recovered Cis is rather poor compared to what should 
ultimately be achievable for a 15°-radius survey. The effect of this is that the recovered Cis accurately follow the true C\s 
only for / > 50. 



2 The criterion for invertibility is equivalent to being able to estimate the correlation functions for all angular separations (Chon et al. 
2004). A sufficient condition is thus that the observed region has pixels separated by all angles between and 180°. 
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3 EXACT COVARIANCE PROPERTIES 



The covariance of any quadratic power spectrum estimates that derive from the pseudo-Cjs by linear transformations follow 
simply from the covariance of the latter. For this reason, in this paper we concentrate on the pseudo-C; covariance. Furthermore, 
we focus on the sample covariance (i.e. we neglect instrument noise) to highlight those aspects specific to CMB polarization 
that arise from mode-mixing on the incomplete sky. Including instrument noise is obviously important for application to near- 
future surveys since it will dominate sample variance. Extending the analytic results derived in this paper to include simple 
white noise should be straightforward (see Efstathiou 2004 for the case of the temperature anisotropies) . However, for more 
realistic correlated noise, Monte-Carlo simulations will probably be required to account properly for the noise contribution 
to the covariances. We shall assume throughout that the CMB signal is Gaussian. If the initial fluctuations were very-nearly 
Gaussian, as expected in simple single-field inflation models (e.g. Bartolo et al. 2004), non-linear effects will introduce only 
a small level of non-Gaussianity into E modes on the scales of interest for primary CMB science (e.g. Mollerach, Harari & 
Matarrese 2004). However, if the primordial tensor-to-scalar ratio r < 0.01, then the dominant power in B modes is expected 
to come from gravitational-lensing conversion of Gaussian E modes into non-Gaussian B modes. In this limit, our assumption 
of Gaussianity will be violated and the expressions derived in this paper will require some modification to account for the 
non-zero connected four-point function. While a full treatment of the pseudo-Ci covariance in the presence of lens-induced B 
modes is beyond the scope of the current paper, we note that the lensing contribution to cov(CP , Cff) has been calculated 
recently for uniform weighting in the flat-sky limit (Smith, Hu & Kaplinghat 2004). 

In this section we calculate the exact sample covariance of Cf and Cf (see also Hansen & Gorski 2003 for a similar 
calculation) and present the results of numerical computations for two types of survey: a full-sky survey with a simple Galactic 
cut, and a 15°-radius circular survey of the type that is optimal for gravitational wave searches in the presence of instrument 
noise at the level achievable with current technology (Jaffe, Kamionkowski & Wang 2000; Lewis et al. 2002). 

For Gaussian CMB fields, the pseudo-multipoles are also Gaussian, so the covariance of their power spectra are given by 

e.g. 



cov(Cf,C^ 



lmE*l m y)\ 2 . 



(2Z + 1)(2Z' + 1) 

If we expand the pseudo-multipoles in terms of the true multipoles using equation (3), we find 



(24) 



cov(c ( ,cv; 



cov(C; ,C V 



(2Z + 1)(2Z' + 1) 



(2Z + 1)(2Z' + 1) 



(2i + i)(2Z' + 1; 



E 

mm' 

E 

mm' 

E 



E +I(l™)(LM) + I(lm.y (LM) 



cE + 



'(lm)(LM) 



E 

LM 

E 



l(lm)(LM) 



-I(lm)'(LM)C'L 



(lm)(LM) 



-I(lm)'(LM) 



-I(lm)' (LM) 



cE 



cE 



+ I(lm)(LM)-I(l m )' (LM)Cl + 



I(lm)(LM) +I(lm)' (LM) CE 



(25) 



(26) 



(27) 



Evaluating these expressions numerically is straightforward in principle, but is computationally expensive, with an operation 
count of 0(l 6 ) if multipoles up to I are retained. This is prohibitive for Z above a few hundred, and so, in general, there is a need 
for faster approximate calculations. One such approximation is to compute the covariances from Monte-Carlo simulations; this 
has the benefit of allowing one easily to include a number of real- world effects that are more difficult to include in the analytic 
calculation. An alternative, which is developed here in Section 4, is to use analytic approximations (that can be calibrated 
with simulations if required). As we show later, it is difficult to find good analytic approximations at low I because of the 
expected sharp rise in the polarization Cis there, and the increased importance of E-B mixing. A similar problem arises for 
the temperature anisotropies because the C;s do not go smoothly to zero as I — > 2. Efstathiou (2004) suggests replacing the 
block of the approximate covariance matrices with small I and I 1 with the exact expression, since the latter can be calculated 
quickly at small I and I' . For diagonally-dominant matrices, this procedure returns an accurate covariance matrix; a similar 
procedure can be implemented for polarization using equations (25)-(27) and the approximations developed in Section 4. 



3.1 Examples 

As our first application of equations (25)-(27), we consider a full-sky survey from which we remove a ±20° band around the 
Galactic plane. We adopt a north-south symmetric weighting which is uniform (which we expect to be nearly optimal in the 
noise-free limit) except for the edge of the retained survey which we apodize with the square of a cosine. In the northern 
hemisphere 

fa a\ - J 1 for 6 < 65 ° e>*\ 

w ^^>-\ cos 2 [36(7O°-0)] for 65° < < 70°. 
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50 100 150 200 

I 

Figure 2. The covariance matrix cov(Cp , Cff) for a Galactic cut of ±20° apodized with the square of a cosine over a width of 5°. The 
oscillatory structure at low / (which saturates the colour scale) is due to the reionization feature in and residual Fourier ringing in 
+I(lm)(lm)' ■ in this plot, as with all plots of covariance matrices in this paper, the units are fj,K 2 . 

Since w(n) is azimuthally-symmetric, the matrices ■izla m )a m ') are block diagonal and the covariances can be evaluated exactly 
in 0(l 4 ) operations. 

Figures 2-4 show the exact covariance matrices cov(C'f', Cfi), 1(1 + 1)1' (I' + l)cov(Cp , Cjf) and W cov(Cp , Cfi) respectively. 
The prefactors have been chosen to make the elements near the diagonal more nearly uniform. Given that the weight function 
is symmetric under parity, +I(j, m )(i m y = +I(im)(i m y and -I(i m )(i m y = -I(i m )(i m y, so that cov(Cf , Gp ) 

and cav{C?,C$) vanish if / + /' is odd, while cov(C , i B , Cy~) vanishes if I + l' is even (and hence on the diagonal). In each of 
the three covariance matrices, the intense feature for I and l' < 50 (which saturates the colour scale) arises predominantly 
from the additional power in Cf on large scales due to reionization. (It vanishes if we substitute for power spectra with no 
reionization.) 

In the case of cov(C l i B , Cjf), the first term inside the summation over L and M in equation (25) is always dominant. 
The covariance is thus very similar to that for the temperature anisotropies with the same weight function. The matrix is 
diagonally dominant, except on large scales, with width inversely related to the angular dimension of the retained region. 
On large scales the covariance matrix is dominated by the additional polarized power generated by reionization. Note that 
this power propagates to somewhat smaller scales in the pseudo-C; covariance due to residual Fourier ringing in the matrix 
+I(im)(im)' ■ (For the weight function adopted here, the width of the main peak of ±I(i m )(im)' is given by the inverse of the 
survey size, but the effective band-limit is larger being set by the inverse size of the apodized region.) 

To emphasise the importance of E-B mixing in the structure of the covariance matrix for Cf , in Fig. 3 we show both 
the exact cov((7 i B , and the matrix we obtain by setting Cf = 0. The latter is very inaccurate on the diagonal at I < 120, 
where E-B mixing is most important, and particularly so at the largest scales where it further fails to pick up the additional 
power in Cf from reionization. Setting Cf to zero also artificially suppresses the correlations: there is significantly more 
structure off the diagonal in the left-hand plot in Fig. 3 than on the right. This reflects the fact that -I{i m ){im)' is intrinsically 
broader than +I{i m )(i m ') (see Fig. 1). 

The covariance matrix cov(C' i E , Cjj), shown in Fig. 4, is only non-zero because of E-B mixing. It is correspondingly 
broader than cov(Cf , Cy). The dominant contribution is from i?-mode power, so the covariance matrix peaks near the 
diagonal at the positions of the acoustic peaks in Cf. 

As our second example we consider the 15°-radius survey with cosine apodization over the last 5° that we used in Fig. 1. 
We should expect the relative importance of E-B mixing to extend to higher I for this smaller survey, and this is indeed seen 
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Figure 3. The covariance matrix /(i + l)/'(/' + l)cov(C ; B , C^) for the same Galactic cut as in Fig. 2. The left panel is the exact covariance 
matrix, while the right panel has Cp = to show the importance of E-B mixing for this covariance matrix. We have multiplied by 
1(1 + + 1) to make the diagonal elements more nearly uniform. The oscillatory structure in the matrix on the left at low I is due 

to the interaction of the reionization feature in Cp with residual Fourier ringing in the coupling matrix —I(i m )(lm)' • 





I 



300 



400 



500 



Figure 4. Blocks of the covariance matrix ii'cov((7 ; B , C^) for the same Galactic cut as in Fig. 2. We have multiplied by IV to make the 
elements near the diagonal more nearly uniform. 

in Fig. 5 where we compare the exact cov(C , ; s , Cjj) with the covariance obtained by setting Cf = 0. The iJ-mode power 
dominates the covariance for all scales plotted (maximum I = 500), and the magnitude of the diagonal follows the acoustic 
peaks in Cf . Surveys of this size are being considered by a number of upcoming experiments targeting B-mode polarization 
from gravitational waves. It is clear that, even for the relatively optimistic tensor-to-scalar ratio adopted here (r = 0.15), 
the application of pseudo-C; techniques to such surveys will demand careful modelling of the impact of E-B mixing in the 
covariance matrix for Cf on any scales where sample variance is expected to dominate. 



4 APPROXIMATE COVARIANCES FOR SMOOTH WEIGHTING 

Given that exact evaluation of the pseudo-C; covariance matrices is, in general, prohibitively expensive, in this section we aim 
to derive accurate analytic approximations that allow fast computation of the covariances at high I. This problem has been 
well-studied for temperature anisotropies (Hinshaw et al. 2003; Chon et al. 2004; Efstathiou 2004), but, as yet, there has been 
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Figure 5. Blocks of the covariance matrix 1(1 + 1)1' (I' + l)cov(C/ 5 , C ; f ) for a 15°-radius region with weighting as in Fig. 1. The exact 
covariance matrix is shown on the left, and the contribution from B modes is shown on the right (i.e. with C ; B set to zero). 

no comparable treatment for polarization. As with the temperature case, the main simplification arises from assuming that 
we are working at sufficiently high I, and with a weight function that is sufficiently band-limited, that we can approximate 
the Cis as being constant over the width of all coupling matrices. This removes the C;s from all convolutions. For +I(i m )(im)' , 
and its products, this approximation should be good, as in the temperature case, for smoothly-apodized observations covering 
a connected patch of linear extent more than several degrees. The approximation is poorer for -Ia m )a m y i an< A breaks down 
completely if there is no apodization, i.e. the weight function w(n) is set to unity or zero. (See Lewis et al. 2002 for examples 
of elements of -In m )a m y m this case.) Here we restrict ourselves to situations where these assumptions hold. We first consider 
a specific case of Gaussian weighting over small patches of the sky, before considering general (but smooth) weighting on the 
sphere in Section 4.2. 

4.1 Gaussian weighting in the flat-sky limit 

As a warm-up to the approximate calculation of the pseudo-C; covariances for general weighting functions w(n), we first 
consider Gaussian weighting over small patches of the sky. For such observations we can approximate the spherical-harmonic 
analysis by Fourier analysis, and we are able to make analytic progress more simply than for the general case. The results of this 
subsection are also of considerable interest for observations of CMB polarization with interferometers [such as DASI (Kovac 
et al. 2002), CBI (Readhead et al. 2004) and the Array for Microwave Background Anisotropy (AMiBA) 3 ], since they directly 
sample the Fourier transform of the sky multiplied by the primary beam. The latter thus acts as a weight function w(n) for 
an interferometer. Furthermore, the primary beam can often be approximated as a Gaussian. 

In the flat-sky approximation, for lines of sight close to the z axis, we adopt a global polarization basis that is aligned 
with the x and —y axes. The radiation propagation direction at the centre of the field is along — z. Stokes parameters defined 
on this basis, Q(x) and U(x), at position x on the sky can be decomposed into Fourier modes Q(l) and U(l): 

Q(l) = j^-Q( x )e- U *, U(l) = J^-U(x)e- il *. (29) 
3 http:/ /amiba.asiaa/sinica.edu.tw/ 
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From these we can define the Fourier transforms of the electric and magnetic parts of the polarization, 

-(ETiB)(l) = e ±2 " t ''(Q±iU)(l), (30) 

where <j>i is the angle between I and the x axis. Note that this is a rotation of the Stokes parameters in Fourier space onto a 
basis adapted to the wavevector I. The Fourier transforms E(l) and B(l) can be shown to be related to the spherical multipoles 
Ei m and Bi m for I >• 1 by the same correspondence as for the temperature anisotropies, i.e. 

E lm « (~i) m \[^ J d4>i E(l)e-' m ^ . (31) 

A similar relation holds for B. In terms of the power spectra, we have e.g. (E(l)E* (V)) = Cf8(l — I'). 

If the polarization is weighted by w(x), then the Fourier transform of the Stokes parameters for the weighted field (or 
the visibilities in the case of an interferometer) are 

0(0 = f ^-w( X )Q( X )e~ il x , (32) 
£L w (l-l')Q(l'), (33) 



/ 



with a similar expression for U(l). Here w(l) is the Fourier transform of w(x). For w(x) real, Q*(l) = Q(—V) and similarly for 
U(l). Decomposing Q(l) and U(l) as in equation (30), we find 

(-E±iB)(l) = f ^-e ±2i{ * l -^' ) w{l-l')(-E±iB)(l'), (34) 

J 27T 

or equivalently 

E(l) = J ^[ + I(l,l')E(l')+i.I(l,l')B(l')}, (35) 

B(l) = j ^[+I(l,l')B(l')-i.I(l,l')E(l')\, (36) 
where the Hermitian [for w(x) real] kernels are 

+1(1,1') = w(Z-Z') cos 2(0, -1(1,1') = -iw(l-l')sm2(<t>i-(j>i'). (37) 

These are the flat-sky limits of the ±Iam)(lm)' m equations (6) and (7). The geometric mixing of e.g. E(l) into B(l) is suppressed 
by the trigonometric factor sin 2(0; — <?V). If the extent of the support of w(l) is ~ / max [i.e. the band-limit of w(n)] then the 
mixing of E into the observable B will be ~ (lmax/l)E(l) for I i ma x- Clearly this is significant if \B(l)\ < (l ma x/l)\E(l)\. 



4-1.1 Pseudo-Ci covariance 

In the flat-sky limit the pseudo-CiS are defined by 



C 



\E(l)\ 2 



C 



-AS 



d<f>i\B(l)Y 



(38) 



The signal correlations between the E(l) and B(l) are discussed in Appendix A. For Gaussian weighting, w(x) = cxp(— x 2 /2a 2 ), 
and smooth spectra (over a range of multipoles ~ l/o"), we can use equations (A1)-(A3) to show that 



(Cf)^a(l)Cf+(3(l)CF, 
whore 



«(0 = 


a 2 
2 


1 1 

2 ~ Pa* 


(3(1) = 


1 

2P 


[l 3 



(C?) ^ a(l)C? + f3(l)CF , 



1 + 



2Z 2 a 5 



2l 2 a 2 



i+ 



Per 2 



(39) 

(40) 
(41) 



These results are consistent with equations (15) and (16), if we use J2i' Pw ~ Q (0 an d J2i> ~ 0(0 f° r ' ^ 1; since the 
(spherical) power spectrum of the Gaussian weight function wi w 7ra 4 cxp(— £ 2 <r 2 ). 

The sample covariances of the pseudo-Czs are given by the angular average of the squares of the correlators of E(l) and 
B(l): 



cov(Cf,Qf) 



- \! 



d^d<f> v \(E(l)E*(l'))Y 



(42) 
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cov(Cf,Cf) 



COv(Cf,Cf) 



dfrdfr,\(B{l)B*{l'))\ 2 , 
dfrdfr,\{E(l)B*(l'))\ 2 . 



(43) 
(44) 



For an azimuthally-symmetric weight function, the correlators in these integrals are invariant under rigid rotations of I and 
I' , so the integration over one of (pi and (j>p is trivial. To evaluate equations (42)-(44) for the case of Gaussian weighting we 
make use of the results for the correlators given in Appendix A. The analysis is rather involved so we simplify things a little 
by considering only the limits Cf = and Cf = 0. If we consider the covariance of Cf with itself, the case Cf = describes 



the contribution to the covariance if there were no E-B mixing, while the case C x 
solely from mixing. 

Consider Cf = first. Then, for azimuthally-symmetric w(x), we have e.g. 



describes the contribution that arises 



cov(Cf,Cf) 



cfcf 



dip 
2^ 



^ + i(i,L) + r(i',L] 



(45) 



where ip = fai — 4>i, and similar expressions hold for the other covariances. We have not been able to perform the integrals 
over ip exactly, but it is possible to make progress in the limit of large I and I' (>• 1/a). In this limit, we can ignore the 
terms with factor exp[— (I 2 + l' 2 )a 2 /2] in equations (A3) and (A6) in Appendix A. This leaves a factor of exp(— \l — l'\ 2 a A /2) 
multiplying slowly-varying functions of ip. The exponential peaks sharply around ip — (reflecting the fact that E(l) and 
B(l') arc only strongly correlated for I within a radius ~ 1/a of I'), so we can accurately handle the inverse powers of \l + 1'\ 
and trigonometric terms in equations (A3) and (A6) with an expansion in sin 2 ip. To make this expansion it is useful to note 
that 

,2 Y 2 )2 



COS 2((pl + <j>i' - 2(pi +l r) 



1 - 2 



sm2(<t>i + fa, - 2<p l+l ,) = 2 



(I 2 



\l + l 

-I' 2 ) 



'14 



■ sin ip, 



\l + V 



[(I 2 + I' 2 ) cos t(j + 211'] sinV>. 



(46) 



(47) 



Integrating the square of the expansions over ip then yields a series of modified Bessel functions hnill'o- 2 ). These series are 
cumbersome, so we only give the leading-order results in the asymptotic regime (I, I' >• 1/a), where we can replace the Bessel 
functions with their asymptotic expansions: 



-(l-l') 2 a 2 /2 



cov(C ! 



8 V2nll'a 2 
6 CfC ( ? 



cov(C, ,C(/ 



{l + l'Y V2nll'a2 
la 2 CFCi 



-(l-l'fcr 2 /2 



(Cf =0). 



(48) 
(49) 
(50) 



2l'(l + l>) 2 s/zrtW 

To obtain these results it is necessary to expand the integrand in equation (45) to 0(sin 4 ip). It is straightforward to verify that 
the covariances of the pseudo-Czs for the other limit, Cf = 0, can be obtained by interchanging E and B in these expressions. 

To gain an indication of how the covariance of the pseudo-C;s propagates to estimates of the power spectra, we form 
simple estimators for Cf and Cf that are local in Cf and Cf by solving equation (39): 

1 



Cf = 



[a(l)Cf - f3(l)Cf], 



Cf = 



[a(l)Cf ~ /3(l)Cf], 



(51) 



iV(I) rr '"' " ~ l N(l) [ 

where the normalisation N(l) = a 2 (l) — P 2 (l). The covariance of these estimators then follows from the covariance of the 
pseudo-C;s which we derived above. To leading order, we find 

6 1 



cov(Cf,Cf 



and 



cov(Cf, Cf 



l&CfCf p _ ( ,_,> ) a g a /a 
V2TTll'a 2 



W + l') 



(l + l') 2 a 4 



I 



I' 



j,3 + ,3 



1 



(Cf=0), 



2CfCf 



-(;-;') 2 <r 2 /2 



(Cf = 0). 



(52) 



(53) 



V2nll'a 2 

Note that these expressions are symmetric in I and I' as required. Equation (53), which gives the covariance of the estimated B- 
mode power if there were no E-B mixing, is the direct analogue of the asymptotic expression for the temperature anisotropics, 
where the exact result in the flat-sky limit (for smooth Cf ) is 

cov(Cf , Cf,) = 2CfCfe-^+ l ' 2 ^^ 2 Io(ll'a 2 ). (54) 

If we average the estimators into bands of width Al >• 1/a we obtain quasi-uncorrelated estimates. The variances of these 
band-powers follows by integrating equations (52) and (53) over I and I' within the band. For wide bands, we find the 
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approximate results 

ftf'E 2 

var(A S ) « (Cf = 0). (56) 

For the Cf' = case, we recover the approximation given for temperature by Hivon et al. (2002), 

var(C, ) « ! (57) 

(2« + l)A// sky w(2) 2 ' 

noting that, for the case of Gaussian weighting considered here, w^'/sky = cr 2 /4 and w^'/sky = cr 2 /8. In Section 4.4 we 
develop more general rules of thumb for arbitrary (but smooth) weight functions; these properly reduce to the expressions 
derived here for Gaussian weighting in the flat-sky limit. 

Our leading-order results show that, for Gaussian weighting, the ratio of the cosmic variance contributions to Cf from 
E and B is ~ (Cf /Cf) 2 /(la) 4 for I ^> 1/a. In Section 5 we consider the implications of this residual variance from £-mode 
power in unbiased quadratic estimates of Cf for detecting gravitational waves via B-mode polarization. 

4.2 Approximations for general weighting 

We now consider the case of general (but smooth) weight functions and work on the spherical sky. The starting point is to 
approximate the correlators of the pseudo-multipoles that appear squared in the pseudo-Cj covariances, equations (25)-(27), 
by removing the C;s from the convolutions: 



(E lm E* lm y) — ^2 +I(lm)(LM)+I*lmy (LM)Cl + -I(lm)(LM) -I(lm)' (LM)Cl 
LM 

« ^CfCf[ + I(w% lm y lm y - [jcfCf - y/cfcfj [-l\u,)] (lmWm y, (58) 

(Bl m B* lm y) = +I(lm)(LM)+I*lm)' (LM)Cl + - I(lm)(LM) - I(lm)> (LM)Cl 

LM 

« JcfCf[ + I(w 2 )] {lm){lm y + ^CfCf - yJcfC^ \-I 2 (w)] {lmWmy , (59) 

{El m B* lm y) = +I(lm)(LM)-I(lm)' (LM)Cl + -I(lm)(LM)+I(lm)'(LM)C'L 

LM 

« \ (^CfCf + ^CfCf^ [-I(w 2 )] ilm)(lmy 

+ \ (sfcfCf - sjcfCf) [+I(w)-I(w) - -I(w) + I(w)] {lm)(lm y , (60) 



where [±I(w 2 )]( lm y lm y are the matrices ±I(i m )(lmy but evaluated with weight function w 2 (h) rather than w(n), and 
[_/ 2 (w)]( im )( im )/ is the matrix product J2lm -I(lm)(LM)-I(LM){im)' - 4 We now consider the derivation of equations (58)-(60) 
in detail. We start from the completeness relation J2i m sYi m (n) s Y* m (h') = S(n — n') which ensures that 

[±2l 2 (w)] (lm)(lm) , = J2 ±2l ( lm )( LM )±2l(l m )'(LM) 



LM 



Ah w 2 (n)±2Y l * m (n)± 2 Y ilrn y(n) 



= [±2l(w )]( lm )(i m y. (61) 

For w(n) = 1 or 0, so that w 2 (n) = w(n), these relations imply that ±2l(i m )(im)' are projection operators (Lewis et al. 2002). 
If we express ± 2 I in terms of ±1, using the obvious matrix notation, we find 

± 2 I(w 2 ) = + I 2 (w) + -I 2 (w) ± +I(w)-I{w) ± -I{w)+I{w). (62) 
Adding and subtracting gives 

+I 2 (w) + -I 2 (w) = +I(w 2 ), +I(w)-I(w) + -I(w)+I(w) = -I(w 2 ), (63) 
4 Due to the Hcrmiticity of ±I(i m )(i m y , [-I 2 (w)] {lm){lm y is also equal to Y.LM - I (lm){LM) - J (*i m )'(LM) • 
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which can be used to verify the second equalities in equations (58)-(60). Note that if Cf were equal to Cf , the correlators 
would only involve ±I(w 2 ) and the construction of the covariances would then be no more difficult than for the temperature 
anisotropies. The major simplification that arises in this case is that the correlators of the pseudo-multipoles reduce to linear 
functionals of w 2 (n). For statistically-isotropic CMB signals, the pseudo-C; covariancc must be invariant under rotations 
of w(n), so that if the correlators are linear in w 2 (n) then the covariance matrix can only depend on rotationally- invariant 
quadratic combinations of the multipoles of w 2 , i.e. the power spectrum of w 2 . In the general case, where the pseudo-multipole 
correlators are more-general quadratic functionals of w(n), the covariance matrix can depend on less-specific configurations of 
the trispectrum of w(ii). (The power spectrum of w 2 (n) picks up only rather specific - and easy to calculate - configurations 
of the full trispectrum.) 

Forming the covariances from equations (25)-(27), we now find 

: + 1) L + 1) E y/c^[+I(w 2 )] {lmKlm y - UC?CB - y/cpcA \-I 2 (w)] {lmWm y ' , 



(<5f,<5 ; ?) 



(2Z + 1)(2J 
2 



(2/ + l)(2i' + l) 
2 



TTjTTJ E yC?C*l + I(w% lmKlm y - (yJC?C» - yJCfC^j [-l\w)] {lm){lm y 



(64) 
(65) 



(2/ + l)(2i' + l) 



\(lm)(lm)' 



+ l - [yJCfCE- y/CfCf) [+I(W)-I( W ) - -I(W) + I(W)] 

It is straightforward to compute ±I(w 2 ) in terms of Wigner-3j symbols using equation (5): 



J(w 2 )] 



(lm)(lm)' 



(2Z + 1)(2Z' + 1)(2L + 1) 



4% 



[1±(-1) K ] 



I I' L 

-2 2 



/ 

m 



(lm)(lm)' 
I' 

-m' 



L 

-M 



(66) 



(67) 



where (w 2 )lm are the (scalar) multipoles of w 2 (n), and, as before, K = I + I' + L. For -I 2 (w), we go back one step, taking 
as our starting point the integral expression for -I(i m )(im)'- 

-I(lm)(lm.y = y dn w(n)| [2YC m (n)2Y { imy(n) - -^^(h) - 2 Y {lm y (n)] . (68) 

Our aim is to find an approximation to _/ 2 (ui) that is a functional of the squares of derivatives of w, so that its contribution 
to the covariance matrices will only depend on easily-calculable configurations of the trispectrum of w. Writing the spin-2 
harmonics in equation (68) in terms of spin-weight derivatives of the scalar harmonics (see e.g. Appendix B of Lewis et al. 
2002), 



iYlm = 



0-2)! 



(1 + 2)! "'" 
which follow from repeated use of 



2Y, m = 



(i + 2)!° ;m ' 



dsY lm = y/(l-s)(l + s + l) 3+ iY lm , d s Y lm = + s)(l - s + l) 3 -iY lm , 

integrating by parts and using (9 2 9 2 — 9 2 9 2 )Yj m = 0, we find that 



(69) 



(70) 



l (lm){l m y 



(l' + 2)\ 



J dn [b s2 w(n) 2 Y* m (h) - b s2 w(h)- 2 Y l * m (h) + 2o 1 w(n)b s 2 Y* m (h) - 2dw(n)Q- 2 Y* m (n)] F (im )/ (n).(71) 



This result shows explicitly that -I(i m )(lm)' vanishes if w(n) is a constant over the full sky. More generally, for w(n) band- 
limited to / m ax, it suggests that -I(im)(lm)' ~ O(Z max /0 + ^( im )( i '»)' f° r ' ^ max with this leading-order contribution coming 
from the last two terms in equation (71). Again, this is consistent with our expectation that E-B mixing is suppressed for 
I ^> ^max- Forming the square of -I(i m )(i m y , we find 



./ 2 (w)]( !m )( im )/ 



i J dnidn 2 | [d 2 w(n 1 ) 2 Y l * m (n 1 ) - b a w(hi)- 2 Y l "" m {h 1 ) + 2Qw(n 1 )B 2 Yi* m (n 1 ) - 2b i w(n 1 )b i - 2 Y{' m (n 1 )] 
x [d 2 w(n 2 ) 2 Y ( * lm y(n 2 ) - d 2 w(n 2 )- 2 Y ( * lrn y(h 2 ) + 2dw(n 2 )d 2 Y ( * lm y (n 2 ) - 2Bw(n 2 )B- 2 Y ( * lm y (n 2 )] * 

(72) 



E Y LM {h,)YZ M {h 2 ) { ±-^\ 



Although we shall not make use of it here, we note for completeness that the sum over L in this equation can be performed 
exactly as follows: (i) apply the operator V 2 (V 2 +2) to remove the factor (L — 2)!/(L+2)!; (ii) perform the resulting summation 
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with the addition theorem; (iii) solve the differential equation away from fi\ ■ fii = 1; and (iv) enforce the correct behaviour 
at h\- fi2 — 1. The final result is 



l-21n2 l + 61n2 „ 1 ,„ mi ,„ 

i67r + 4 8 ^ co S/ 3+— (1- cos/3) ln(l- cos/3), 



(73) 



where cos/3 = fii ■ n 2 . The first two terms can be shown not to contribute to -I 2 (w). For our purposes, it is more useful to 
approximate the sum in equation (72), for 1,1' ^> i m ax, by 



E (L + 2)! YLM ^) Y LM{n2) 



(1-2)1 (l'-2)\ 



(l + 2)\V (l' + 2)\ 



E Y LM (fii)Yl M (n 2 ) 



(1-2)1 (l'-2)\ 



(l + 2)\V («' + 2)! 



5(ni 



n 2 ) 



(74) 



Note that this expression only makes sense when considered inside the integrals in equation (72). The justification for the 
approximation is as follows. Equation (72) may be regarded as a convolution of the slowly-varying (L — 2)\/(L + 2)! with the 
product of two functions of L: one resulting from the integral over ni, which is localised to within i max of I, and the other from 
the integral over n 2 , which is within £ max of I'. We can thus replace (L-2)\/(L + 2)\ with yf(l - 2)\/(l + W-VW - 2 )-/( 1 ' + 2 ) ! 
in the convolution for I and I' S> ^max- Furthermore, we can then extend the summation to include the L = and L = 1 terms 
since they sum to a constant plus cos/3 which, like the first two terms on the right of equation (73), do not contribute to 
the integral. (These terms would be negligible anyway, for large I and I' .) The presence of the delta function in equation (74) 
ensures that we achieve our aim of reducing _/ 2 (to) to a linear functional of squares of derivatives of w(n). However, before 
proceeding we make one further approximation which can easily be relaxed if more accuracy is required. Since w(n) is assumed 
to be slowly varying compared to the spherical harmonics, we ignore terms like d 2 W2Y l * m compared to e.g. QwB 2 Y l '$ H , so that 
we retain only the leading-order contribution to -I 2 (w). Finally, with these approximations we have 



-I 2 ( w )](lm)(lm)> 



(/ + 2)i y (I' + 2)! / d " {[^ w ( A )^ Y *m(n) - 3™(n)5L 2 Yi*m(n; 
x [dw(n)d 2 Y( lm y(h) - fiw(n)B-2Y( lm y(h)]*} 

(_ 1 )m+l 



y/l(l + l)l'(l' + l) 



E 



(2/ + i)(2Z' + l)(2L+i) [ I I' 

rn —m' 



x |3< M [1 + (-!)*] 



I V 
-1 1 



47T 

L 




L 
-M 



+ [ 2 (9w)!m + (-l) f 



-2(9w)| M 



-1 



(75) 



where we have used equation (70). Note that this approximation preserves Hermiticity of -I 2 (w). In deriving equation (75) 
we have expanded products of the derivatives of w(n) in the appropriate spin-weight harmonics: 



Qw(n)[dw(n)]* = dw(n)[dw(n)]* = \Qw(n)\ 2 = E \^ w \tmYim(n) , 
dw(n)[Bw(h)]* = [dw(n)] 2 = ^ 2 {3w) 2 m2 Yi m (n) , 

1^2, m 

Qw(n)[dw(n)]* = [dw(n)] 2 = ^ -2{Sw) 2 m - 2 Yi m (n). 

1^2, m 



(76) 
(77) 
(78) 



Note that (5w(n) is the complex conjugate of <5w(n) since w is real, and that the spin-2 multipoles satisfy _2(3w)? m = 
(— l) m 2(9w)f* m . Substituting equations (67) and (75) into equations (64) and (65), we have 



cov(cf,Q?) = -L^J [! + (-!)*] /^(» 2 ) 

LM [ 



y/l(l + l)l'(l' + l) 



I I' 

-2 2 

I V 
-1 1 



+ Elm 



^/l(l + l)l'(l> + 1) 



E[i-(-i) K ]|eLM| 



(79) 
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and 



LM I 



CfC^(w 2 ) LM 



+■ 



y/l(l + l)l'(l' + l) 

^Z(Z + 1)Z'(Z' + 1) 



3w|lm 



Z Z' L 

-2 2 

Z Z' L 

-1 1 



Z Z' L 
-1 -1 2 



£[l-(-l)*]|£ iM | 5 



Z Z' L 
-1 -1 2 



(80) 



where we have introduced parity eigenstates Ei m and Bi m for the multipoles of the spin-±2 objects [3to(n)] 2 and [dw(n)] 2 , 
which are defined by 5 

2 (9w)? m = Elm + iBim, -2(5w) 2 m = £l m - iBl m . (81) 

Note that the Bi m vanish for an azimuthally-symmetric weight function since then 2(3w)?o = -2(8w)f . Equations (79) and (80) 
achieve our goal of finding approximations to the pseudo-C; covariance matrices that are efficient to compute, but do take 
account of E-B mixing effects at leading order for I, I' 3> Z max . Note that, whereas the pseudo-C; covariance for temperature 
anisotropies involves only the power spectrum of w 2 (e.g. Efstathiou 2004), the polarization covariances also include the auto 
and cross-power (with w 2 ) of the square of the gradient of w. 

We can make a useful check on our approximation for {-I 2 (w)](i m )(im)' , equation (75), by setting Z = Z' and m — m! and 
summing on m. Inspection of equation (16) shows that 



ZjT[ I][- j2 MW)(;m) = M »'' 



2Z + 

ui, i ■ 

whereas evaluating the left-hand side by summing the approximation for _/ 2 (to) directly gives 



2Z 



V^Ki + 1) 



(|3H ; 



\X7w\ 2 
2tt J 1(1 + 1 



■ dn. 



(82) 



(83) 



Reassuringly, this agrees with the leading-order result for J2i> given in equation (19). 

Finally we consider the covariance of Cf and Cj? , equation (66). This involves evaluating the matrix commutator [+/, _/]. 
Following the procedure that led to equation (71) for -I(i m )(lm)' , we can express +I(i m )(lm)' as 



+ -l(lm)(lm)' — 



(I' - 2) ! 
(Z' + 2) 



j J dn | i [d 2 w(n) 2 Yi m (n) + d 2 w(n)- 2 Yi m (n)] + [dw(h)'d 2 Y lm (n) + Qw(n)Q- 2 Y lm (r 



(1+31 

{1-2)1 



w(n)Y lm (n) > Y (lm y(n 



(84) 



This is similar to equation (71) except for the presence of the last term in the curly brackets (and the sign differences). The 
last term will dominate the integral for Z, Z' 3> Z max , and in this limit +I(i m ){im)' will be close to the equivalent coupling matrix 
for the temperature anisotropies (see e.g. Hivon et al. 2002). Forming the matrix product +I(w)-I(w) from equations (71) 
and (84), and using the approximation in equation (74), we find at leading order 



(2Z + 1)(2Z' + 1)(2L + 1) 

47T 



1 / L(L + 1) 2x i . . ' ' ! - 

2\ Y(FTi) {w)LM[1 - { - l) 1 o i -i 



rn -m' -M 



I V 



(85) 



which is 0(Z max /Z) as expected. In deriving equation (85) we have used (w3w) = 3iu 2 /2 = J2lm V L(L + T)(w 2 )lmiYlm/2. 
Noting that -I(w)+I(w) is the Hermitian conjugate of + I(w)-I(w), we have 



./(.mm - (-i)-Ev/ (2,+i)(a, 4 : i,(2t+i) (,: _'„. 

LM \ 
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-M 



2]j 1(1 + 1) 

5 Under a parity transformation w(h) — » w(—n), 2 (Bw)f m — ► (— I) 1 - 2 (Sw)f m and ~ 2 (3w)f m — » (— l) 1 2 (<5w)'f n 
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(86) 
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where we have used the symmetries of the 3j symbols and the reality of w 2 . Differencing the previous two equations gives 
+I(w)-I(w) — ~I(w)+I(w) 



E(- 



(2J + 1)(2J' + 1)(2L + 1) 



4ir 



w )lm 



x[l-(-l) J 



This can be simplified further by noting that 



[l-(-l) 



-1 



I' 
1 



-1 



1(1 + 1) 



(87) 



,(88) 



which follows from a recursion relation for the 3j symbols [equation (4) of Section 8.6.2 of Varshalovich et al. (1988)]. Since 
\l — l'\ ^ imax, it follows that the commutator [ + I(w), -I(w)](i m )(i m y is O(£ m ax/0 2 > an< A so we can safely neglect its contribution 
compared to that from -I(w 2 ) in equation (66) at high I. (The power spectrum prefactors are approximately equal for both 
terms.) In this limit, cov(Cf , Cf) is symmetric. Therefore our final approximation for the covariance of E and B is, after 
summing over m and m', 

\ 2 

I I' L 



cov (cf ,c?) - ^ El 1 - (-!)"] (7 W + 



IK 



-2 2 



(89) 



This, and our approximations for cov(Cf , Cf ) and cov(Cf ,Cf), can be computed efficiently. For a general weight function 
w(n), we can compute the multipoles |9w| 2 m , £i m and Z3; m as follows. First compute the spin-1 field, 3w; this can be done, 
for example, by synthesising J2i m \/^(^+l) w imiYi m (n) with fast spherical transforms (on an iso-latitude pixelisation) . Next, 
form the spin-zero field |3u>| 2 and the spin-2 field (3w) 2 in real space by taking the squared modulus and the square of 3w 
respectively. The scalar multipoles |3ii>| 2 m can then be extracted from the spin-zero field with a (scalar) spherical transform, 
while £i m and Bi m are just the electric and magnetic multipoles of the spin-2 field which can be extracted with spin-2 spherical 
transforms. At each I and V , the 3j symbols for all the L values required can be computed efficiently with stable recursion 
relations (Schulten & Gordon 1976). 



4.3 Examples 

In this subsection we test the approximate covariance matrices developed above against the exact matrices. Since computing 
the latter on intermediate and small scales is prohibitively slow for general weight functions, we restrict ourselves to the 
two azimuthally-symmetric examples considered in Section 3.1. Of course, for any practical application, our approximate 
covariance matrices should be carefully verified against those obtained from high-volume Monte-Carlo simulations with the 
chosen survey geometry. 

Figure 6 shows the covariance matrix cov(Cf , Cf) for the 15°-radius survey computed with the approximation (79) and 
with the exact expression, equation (25) . Again, we take the weighting to be uniform except for cosine-squared apodization of 
the outer 5° annulus. The approximation is very accurate for the covariance of Cf , except on the largest scales, as expected 
given that our largest potential source of error - the perturbative treatment of E-B mixing - is irrelevant in this case since 
Cf > Cf. For I > 55, the error on the diagonal is better than 2% except around the acoustic trough in Cf, at I ~ 200, where 
it is around —8%. Quite generally, the error on the diagonal is largest at the acoustic peaks and troughs; the approximation 
over-estimates at the peaks, where the curvature of Cf is negative, under-estimates at the troughs, and is best in between 
where the curvature vanishes. This is consistent with the errors that we would expect from removing Cf from the convolution 
in equation (25). The magnitude of the correlations for Cf fall to lower than 1% for \l — l'\ > 60. For I and I' > 55, the 
maximum error in those off-diagonal elements with correlations greater than 1% is ~ 30%, but the error is localised to a 
few elements at the extremes of the diagonal band around the acoustic trough at I ~ 200. This is, however, a very stringent 
error criterion: to get a 1%-correlated element accurate to 30% with Monte-Carlo methods would require ~ 10 5 simulations. 
A more relevant measure for the off-diagonal elements is the absolute error in the correlations, since large fractional errors in 
elements with small correlation are relatively harmless. With this measure, the maximum error in the correlations for I and 
I' > 55 is 0.03. We have made similar comparisons for the case of the ±20° Galactic cut described in Section 3.1. In this case, 
the error on the diagonal is better than 1% for I ^ 55. The only off-diagonal elements with correlations ^ 1% are for I' = Z±2 
and I' = I ± 4, and for all of these the fractional error is less than 2%. 

The same comparisons are made for cov(Cf , Cf) for the 15°-radius survey in Fig. 7. The dominant contribution to the 
covariance matrix over the full range plotted is now from i5-mode power that is mixed into Bi m ; see Fig. 5. As expected, 
the approximate covariance matrix for Cf is less accurate than for Cf since any errors in our modelling of the mixing are 
important for the former except at very large I. The implied improvement in accuracy with increasing I is apparent in the 
figure. More quantitatively, for I > 150 the error on the diagonal is better than 10%. Interestingly, the error follows the same 
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Figure 6. Blocks of the covariance matrix cov(C i jE , Cjf) for a 15° radius region with the same weighting as in Fig 
matrix is shown on the left, and its approximation, equation (79), on the right. 



1 . The exact covariance 



trend as for Cf , i.e. it is largest in magnitude at the acoustic peaks and troughs of Cf, and vanishes at the inflection points. 
This suggests that the main source of error here is from removing Cf from the convolution in equation (26) rather than 
errors in approximating -I 2 (w) by equation (75). (Since -Ia m )a m y is broader than +In m )(lm)' > the fractional error made in 
approximating the convolution is worse for the covariance of Cf than for Cf .) Regarding the off-diagonal elements, for I and 
I' > 70 the error in the correlation is < 0.05 for all elements. The correlation matrix is thus accurately determined except on 
the largest scales, and so, in a practical application, could be used as a template for the shape of the covariance matrix with 
the variances calibrated off simulations. This procedure is similar to that adopted by WMAP for the analysis of the one-year 
temperature data (Hinshaw et al. 2003). The Al = \l — l'\ over which correlations are significant is scale-dependent for the 
covariance of Cf unlike that for Cf; see Fig. 9. The correlations are broadest at the peaks of Cf, where the contribution of 
E modes to the Cf covariance is greatest. This behaviour arises from the different widths of the ±I{im)(lm)' matrices, and 
would be lost in any approximation that did not model the E-B mixing accurately. Repeating the comparison for the case 
of the Galactic cut, we find that the diagonal elements are accurate to better than 10% for I > 50, and better than 1% for 
/ > 150. For the off-diagonal elements, the correlations are accurate to better than 0.03 for I and l' > 55. 

Finally, we compare the exact and approximate cov(Cf , Cjf ) for the 15°-radius survey in Fig. 8. The error on the diagonal 
is better than 10% for all I > 55, and, as above, arises mainly from removing Cf (which dominates) from the convolution in 
equation (27), rather than our ignoring the matrix commutator [+I(w), -I(w)] in equation (66). For parameter estimation, 
the error in the correlation matrix, 

corr(Cf , CfO = cov(Cf , 6$ )/yJvnr(Cf )var(Cf ), (90) 

is the more relevant quantity. The exact correlation matrix is shown in Fig. 9; there are strong correlations on large scales 
due to the power in Cf from reionization, but for / and I' > 55, the maximum correlation is only 0.15 (and occurs on the 
diagonal). Forming the approximate covariance matrix with equations (79), (80) and (89), we find that the error is better 
than 0.03 for all I and l' > 55. For the case of the Galactic cut, the correlations are < 0.01 for I and I' > 55, and the maximum 
error in our approximate correlation matrix is ~ 0.001 over the same range. 

To summarise, the dominant error in our approximate covariance matrices, for the weight functions tested here, is from 
the treatment of the power spectra in the convolutions in equations (25)-(27). The biggest errors on large and intermediate 
scales arise from the terms with the dominant Cf convolved with the broad -Ia m )(lm)' % an< i so are most critical for the 
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Figure 8. As Fig. 6 but for the covariance matrix ii'cov(C' ; B , Cff). 
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Figure 9. The exact correlation matrices corr(C ; E , Cf) (left), corr(C; , Cjf) (middle) and corr(Cf ,Cf) (right) for the 15°-radius 
survey. The colour scale has been chosen so that it is saturated by elements with correlation > 0.35 to emphasise the broadening of the 
correlations for Cf around the acoustic peaks of C, . 

variance of Cf . In forming the correlation matrices, the convolution error largely cancels and so our approximate correlation 
matrices are more accurate than the covariances. For small surveys it may be necessary to re-calibrate our approximate 
variances off simulations. Accurate covariance matrices should then be achievable by combining our approximate correlations 
matrices with the calibrated variances. 



4.4 Approximations for band-power variances 

In this subsection we use our approximate expressions for the covariances of the pseudo-CiS, equations (79), (80) and (89) to 
derive useful 'rules of thumb' for the variances of quasi-uncorrelated band-power estimates of the power spectra. Our results 
generalise equations (55) and (56), derived earlier for Gaussian weighting in the flat-sky limit, and extend the temperature 
result of Hivon et al. (2002), given here as equation (57), to polarization. 

Following our treatment of Gaussian weighting in Section 4.1, we approximate the inversion from pseudo-CiS to estimated 
Cis by the local transformation 

Cf = -^(aiC? - frCf), Cf = ^r(aiG, B - piCf), (91) 

where we have defined ai = ^2 H , Pu> and /3; = M n i. Recall that e.g. ^2 a i Pu/Cf is the contribution of E modes to the 
mean of Cf , and M u /Cf is the contribution from B modes. The normalisation Ni = af — (3? . The covariance of the Cis 
is thus a scaled version of that for the pseudo-CiS. As in the flat-sky limit, we consider only the two extreme cases, Cf = 
and Cf = 0, to keep the algebra manageable. Such limits are sufficient to get the dominant contribution to the variance of 
Cf , and to assess the importance of E-B mixing for the variance of Cf . We only need perform the calculation for Cf = 
since the other case then follows by symmetry. 

We obtain quasi-uncorrelated estimates by averaging the recovered C;s in bands of width Ai > i max , so we begin by 
averaging the pseudo-C; covariances in such bands. Consider first the variance of the binned Cf for I 3> imax- The leading- 
order contribution is from the term involving (w 2 )lm in equation (79). To perform the summation over I' in a band centred 
on I, we first insert a factor (21' + 1)/(2Z + 1), which is effectively unity for l' in the band provided that / ^> Ai. So long as 
we further choose bands that are broad compared to imax, we can reduce the lower limit of the summation over I' to 2 and 
increase the upper limit to oo. (Recall that i max sets the range over which the pseudo-C;s are significantly correlated.) Finally, 
if we approximate Cy as constant over a band we can replace it by Cf , and the summation over the square of the 3j symbol 
can be done in the same way as in the derivation of equations (15) and (16). We find that, at leading-order, 

1 + ff cov(Cf ,Cjf) - -J-^gL£ |(„W = -A- vF>f*,C?C?. (92) 

l'=l-Al/2 LM + 

To compute the variance of the binned Cf, we further average over I in the same band, treating Cf as constant, to find 
var(Cf ) (2r+ 2 1)A ^ (4) / sky Cf 2 (Cf=0). (93) 

This result has the same form as that for the temperature anisotropies. 

For the variance of the binned Cf we must retain all terms in equation (80) that survive setting Cf = 0. We perform the 
summation over V in a similar manner to that described above, making use of some additional results for summing products 
of 3jr symbols that are given in Appendix B. Further averaging / within the band, we find, at leading order, 
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™(Cf ) « (2 *„ + 1)2A ^ E (2 I«m| 2 + \Slm\ 2 + \B LM f) (CP = 0). (94) 



(2/ + i)/2(/+i)2 A r .. , 

This can be simplified further by noting that 

EH 3 H£m| 2 = / |9 W j 4 dn= f(X7w) 4 dh, (95) 
and, from equation (81) and the completeness of the spin-weight harmonics, 

^(\£lm\ 2 + \B lm \ 2 ) = l^2[\2(Bw) 2 LM \ 2 + |-2(8w)Im| 2 ] = / |S»| 4 dn. (96) 

LM LM J 

Putting these results together, we find 

^> w ( 2r + i)Pg + i)^ ^/^ dft (c ^ 0) - (97) 

Note that in the limit C; S = 0, the variance of Cf arises solely from E-B mixing which explains why only derivatives of to 
appear in this equation. 

The final term we require is cov(Cf , Cf). Setting Cf = in equation (89) and using similar approximations to those 
adopted above, the leading-order result is 

cov«7f,<7f) « (2/ + 1) ^ + 1)A ^ E^ + l)K-W 



LM 



i cf 



■ J (Vto 2 ) 2 dn. (98) 

As expected, the relative sizes of var(Cf), cov(Cf , Cf) and vax(Cf) for Cf = are roughly 1, OQ^^/T 2 ) and O^^/l 4 ) 
respectively. 

We now have all the results required to compute the leading-order variances of the band-power estimates for Cf and Cf . 
Combining the (co) variances calculated above with equation (91), we find 

-r(Cf) - , 9 .Jr A> 7 (2)2 (Cf=0), (99) 

~ 2Cf 2 6 / (v^)' 4 ' 2 TO W(V W ) (2)2 4^^. 

var ( C r) ~ /nr I iu» g7? I t \ 9. ^T2-+q o ^73 (Ci =0). (100) 



(2! + l)AI/ sky I 2 (I+l) 2 ^ rf) 2 3 w (2) 4 3 w (2) 

Here, we have introduced the convenient notation AnX^ / s k y = / X 1 dn, where X = w, Vui or idVib in the above, and 
approximated a t w to (2) /sky and A w 2(Vto) (2) / sky /[/(/ + 1)] [see equations (15) and (19)]. 

To obtain the equivalent results for Cf = 0, one should swap E and B everywhere in equations (99) and (100). It is 
straightforward to verify that for a Gaussian weight we recover the variances given in Section 4.1 in the flat-sky limit. We 
now see, quite generally, that the ratio of the excess sample variance from E modes in band-power estimates of Cf to the 
variance from B modes is ~ (C , ! E /C i B ) 2 0(i 4 aax /i 4 ). We consider the implications of this for _B-mode detection with pseudo-Ci 
techniques in the next section. 



5 IMPLICATIONS FOR B-MODE DETECTION 

The realisation that, in linear theory, scalar perturbations do not generate _B-mode polarization (Kamionkowski et al. 1997; 
Seljak & Zaldarriaga 1997) has paved the way for a future generation of instruments designed for high-sensitivity searches for 
gravitational waves via B-mode polarization. Instrument noise can be reduced to the ~ 4/iK-arcmin level by surveying e.g. a 
few-hundred square degrees for around a year with arrays of hundreds of detectors. With instrument noise at this level, the 
sample variance from the B modes produced by gravitational lensing of the scalar E modes dominates the thermal noise in 
the random error budget for Cf , assuming perfect isolation of B modes. (We are assuming the tensor-to-scalar ratio r <JC 0.01 
here, so that sample variance of the linear B modes is sub-dominant.) If we make no attempt to reconstruct and subtract 
the lensing signal, its variance sets a fundamental limit to the one-sigma error of Ar = 1.4 x 10~° for a full-sky survey in 
the null hypothesis, r = 0. This should be compared with the limit of only 0.02 obtainable by analysing temperature and 
electric polarization only. The limit on r is much poorer in this case because of the large sample variance of the dominant 
scalar perturbations. 

If pseudo-C; techniques are to be applied to future, high-sensitivity searches for gravitational waves, it is clearly desirable 
that the imperfect isolation of B modes inherent in such methods should not compromise sensitivity to r. In Fig. 10 we 
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Figure 10. Sample variance errors on the recovered Cf using the estimator of Chon et al. (2004); see also Section 2.1. The survey 
area and weight function is the same as in Fig. 1. The error boxes are the contribution purely from Cf to the one-sigma errors on 
flat band-power estimates with a AZ = 70. These are thus representative of the errors that would be obtained if E and B modes were 
separated (without loss) at the level of the map. They agree well with the simple rule of thumb in equation (99), plotted as the solid 
line. The error bars are the contribution purely from Cp. They agree reasonably with the dashed line, which is the rule of thumb in 
equation (100). Critically, these dominate the errors due to Cf in the angular range relevant to gravity-wave searches with B-mode 
polarization. All errors are computed in the null hypothesis r = 0, so the contribution from Cf arises only from sample variance of the 
lens-induced B modes. 



plot the errors on flat band-power estimates of Cf (with Al = 70) obtained by using the unbiased estimator of Chon et al. 
(2004); see also the discussion in Section 2.1. The weight function applied to the 15°-radius survey, and apodization of the 
correlation function, are the same as in Fig. 1. The errors are obtained from the exact covariance matrices for the pseudo-Cis 
by transforming (two-sidedly) with the appropriate pseudo-inverses P^, 1 and M,7/ . These pseudo-inverse matrices are shown 
in Fig. 1. The errors are calculated for r = in two ways. First, we set Cf — so that we obtain the error in the recovered 
Cf due to sample variance of the lensing B modes. These errors are representative of what would be achieved if we were able 
to separate E and B modes perfectly, and without loss, on all scales in the maps; they agree very well with the rule of thumb 
in equation (99) if we swap the roles of E and B there. We also calculate the band-power errors obtained by setting Cf = 0, 
i.e. removing the B-mode power from lensing; they agree reasonably with the rule of thumb in equation (100). These errors 
arise because, in any realisation, our chosen estimator Cp receives a contribution from both E and B modes, despite only 
coupling to B-mode power in the mean. Such errors can only be removed by isolating B modes at the level of the map, or by 
adopting more optimal, non-local weighting of the data, as in maximum-likelihood power spectrum estimation. Critically, we 
see from Fig. 10 that the errors due to the sample variance of the E modes that leak into the estimator dominate the sample 
variance from lensing B modes on large scales, where the gravitational-wave signal resides. We conclude from the figure that 
pseudo-C; estimators will be highly non-optimal for future experiments with thermal sensitivity better than ~ 4 /iK-arcmin, 
targeting B-mode polarization over a few-hundred square degrees. 6 

We can make a quantitative estimate of the impact of the additional sample variance due to leaked E modes on the 
detectability of gravitational waves as follows. To preserve all information in the pseudo-C;s, we bypass the step of forming 
unbiased estimates of Cf [which involves the singular transformation in equation (22)], and instead estimate r directly from 



° If gravitational lensing were absent, the optimal survey for detecting B-modc polarization with a power spectrum of known shape but 
(low) unknown amplitude is a circle of radius around 10 degrees (Jaffe, Kamionkowski & Wang 2000; Lewis et al. 2002). This assumes 
that the sample variance of the B modes from gravitational waves does not dominate the error budget. If the lens-induced B modes are 
treated as an additional (Gaussian) noise term, the optimal survey area is increased since sample variance becomes an issue again. The 
best compromise is roughly to balance the sample variance of the lens-induced B modes with the thermal-noise variance. As a specific 
example, for 200 pairs of polarization-sensitive bolometers, with each bolometer having sensitivity 200 (iK^/s, and one year of integration, 
a survey covering two per cent of the sky balances the thermal and sample variance. 
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the pseudo-C;s. We concatenate Cf and Cf into a single data vector C\, whose covariance matrix cov(C";, CV) then consists 
of blocks formed from cov(Cf , Cf), cov(Cf , Cf) and cav(Cf , Cf). We assume that all parameters except r are fixed, and 
calculate its one-sigma Fisher error from 

1 ^ d(Ci) d(Cy) - 

(A^ = E dr dr cov (^'). ( 101 ) 

In practice, we invert the covariance matrix with a singular-value decomposition which projects out any linear dependencies 
in the C 1 ; due to the finite sky coverage. For the 15°-radius survey, the one-sigma error on r in the null hypothesis (r = 0) is 
Ar = 0.053. Our sensitivity to r improves dramatically to Ar = 1.096 x 10~ 4 if we artificially set Cf = (and dCf /dr = 0) 
in the analysis. This latter case is representative of what could be achieved with lossless separation of the B modes, followed 
by a pseudo-C; analysis of the Stokes maps of these B modes. In practice, E-B separation is lossy, particularly on the largest 
scales where there is significant _B-mode power from reionization (Lewis et al. 2002; Bunn et al. 2003; Lewis 2003). If we 
assume, crudely, that large-scale modes with I < 20 are lost to the separation, we find a Ar ~2x 10~ 3 for the 15°-radius 
survey. The effect of not isolating the B modes properly in the pseudo-C; analysis thus degrades the limit on r by a factor 
~ 25. 



6 CONCLUSION 

Accurate calculation of the covariance matrix is an important part of power spectrum estimation during the analysis of CMB 
data. A 'correct' estimator (i.e. unbiased), but with poorly quantified errors, is useless for subsequent parameter estimation. 
In this paper we have constructed the contribution from sample variance to the covariance matrices of heuristically- weighted 
quadratic estimates of the CMB polarization power spectra. Exact calculation of these matrices for general weighting is 
numerically feasible at large scales (I less than a hundred or so), but the 0(1 ) scaling is prohibitive on small scales. To resolve 
this problem we developed analytic approximations to the polarization covariance matrices, including the effects of mode 
coupling and, in a perturbative manner, the mixing of E and B modes that inevitably arises with heuristic weighting. Our 
approximations, which generalise earlier work on the temperature anisotropies (Hinshaw et al. 2003; Efstathiou 2004), can 
be computed efficiently to high I. Accurate covariance matrices can be obtained on all scales by replacing the low-Z sectors of 
our approximations with a numerical calculation. 

We tested the approximations against numerical results for two toy-model surveys for which azimuthal symmetry allows 
efficient evaluation of the exact covariances. For the more challenging survey, a 15°-radius survey with cosine-squared apodiza- 
tion of the outer 5°, the most significant errors were on the variance of Cf. These errors arise from our local treatment of the 
i5-mode power spectrum when it appears in a convolution with the broad kernel that describes the mixing of E and B modes 
in equation (26). For the variance of Cf (and the temperature anisotropies) the approximation is more accurate because the 
dominant convolution involves a kernel that is intrinsically more narrow. It will probably be very difficult to improve the 
error on the variance of Cf in an analytic treatment. Fortunately, convolution errors tend to cancel in the approximation 
of the correlation matrices. We observed maximum errors < 0.05 in the correlations, except on the largest scales, suggesting 
that, in a practical application, it should be possible to construct accurate covariance matrices by combining the approximate 
correlations with variances calibrated off simulations. By binning our approximate covariances, we also derived useful rules 
of thumb for the variances of band-power estimates of Cf and Cf . These are particularly simple to evaluate, involving only 
integrals of products of the weight function and its derivatives. 

A further aim of this paper was to quantify the conditions under which pseudo-Cj methods can be applied to gravitational- 
wave searches via B-mode polarization. Measurements of the B-mode power spectrum from gravitational waves with future 
surveys with sub-4 /iK-arcmin thermal noise are potentially limited by the sample variance of lens-induced B modes. We 
showed that for surveys covering one or two per cent of the sky, the additional variance due to leaked E modes dominates 
the lensing noise. Adopting pseudo-Cj methods would therefore by highly sub-optimal: the error on the tensor-to-scalar ratio 
would be ~ 25 times higher than if B modes had been isolated before forming quadratic estimates. Ultimately, this additional 
variance would limit the tensor-to-scalar ratio that is detectable with pseudo-C; methods to ~ 0.05. Map-level E-B separation, 
or more optimal weighting such as in maximum-likelihood power spectrum estimation, will therefore be essential to ensure 
the maximum returns from future B-mode surveys. 
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APPENDIX A: CORRELATORS FOR GAUSSIAN- WEIGHTED E(l) AND B(l) 



In this appendix we derive approximate analytic expressions for the correlations between the flat-sky Fourier transforms E(l) 
and B(l) for a Gaussian weight function w(x). Our expressions are valid for I and I' large compared to 1/a. In this limit, the 
power spectra can be approximated as being constant over the support of w(l — L); see Section 4. Using equations (35) and 
(36) we can approximate the correlators by 



(E(l)E*(f)) 



^c B c* 

2tt 



dJL + I(l,L)+r(l',L) + 



\j c i c v f d 2 L 



2tt 



^-I(l,L).I'(l',L), 



C B C B 2 

(B(l)B*(t)) « V L '' [ ^+I(l,L)+I*(l',L) + v o '_ ' / ^-I(l,L).P(l',L). 
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For a Gaussian weight function we can evaluate the integrals in these equations analytically; we find 
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Note that 

j2 



J %± [+I(l,L)+I\l\L) + .I(l,L).I*(l',L)] = ia 2 cos2(^ - fr,^ 1 - 1 ' ^ /4 , 



.(A3) 



(A4) 



and the right-hand side is +1(1, 1') constructed from w 2 (x). This result follows quite generally from equation (37). To calculate 
the covariance of the errors of the power spectrum estimates we also need the following correlator 



(E(l)B*(l')) 



2ty 



i(i,L)+r(i',L). 
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The integrals in this expression evaluate to 
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We now have 
d 2 L 



I 



— [+1(1, L).P(l', L) + .1(1, L)+P(l', L)] = ~ l -a 2 sin2(0 ; - ^e^-''^ ' /4 , 



(A6) 



(A7) 



and the right-hand side is -1(1, 1') constructed from w 2 (x). Note that (E(l)B* (I')) — if I and I' are parallel or perpendicular. 



APPENDIX B: USEFUL SUMS OF PRODUCTS OF 3J SYMBOLS 

In Sections 2 and 4.4 we make use of a number of results for summing products of 3j symbols. Some of these follow trivially from 
the orthogonality properties of the 3js, but others require more effort. Here, we sketch our derivation of these latter results. 
We make extensive use of well-known properties of the 3j symbols, such as their symmetries and integral representations 
(e.g. Varshalovich et al. 1988). 
In Section 2 we require 

SlS £ ( -i)- (2 r + i)( l _ 2 2 o ) • (Bl) 

It is the presence of the factor of (— 1) K , where K = l + l' + L, that complicates the evaluation of this summation. To proceed, 
we express the product of 3j symbols as an integral of reduced Wigner matrices using 

| 1 i < ni (^ 2n2 (/3)^ 3n3 (/3)dcos/3 = 2^ 1 i ^ ™ 3 ) ( £ £ ^) (for = ru = ), (B2) 

so that 

Si = \ X)(2I' + 1) / d l _ 22 (p)4_2(WooW) dcos/3. (B3) 
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The factor has been absorbed in reducing the product of 3j symbols to a square. The summation over V can be 

performed by first expressing d l 2 _ 2 in terms of d l 22 using the recursion relations for the reduced Wigner matrices (Varshalovich 
et al. 1988), and then using the completeness relation. Full details can be found in Chon et al. (2004), but the final result is 

i ]T(2Z' + 1)4-2 = <5(cos/3 - 1) + csc 2 (/3/2). (B4) 
v 

Inserting this into equation (B3) gives 

Si= ycsc 2 (/3/2)4- 2 (/3)rfoo(/3)dcos/3. (B5) 

Using the differential representation of the reduced Wigner matrices (see Section 4.3.2 of Varshalovich et al. 1988), it can be 
shown that csc 2 (/3/2)d2_ 2 (/3) is a polynomial in cos ft of degree Z — 1. It can thus be expanded as the series Yln=o an( ^oo(P)j 1 - e - 
an expansion in Legendre polynomials. Using the orthogonality of the d l 00 , we see that the integral in equation (B5) vanishes 
if L ^ Z. It can be performed for L < I by replacing d\_ 2 by its differential representation and repeatedly integrating by parts. 
Only boundary terms then survive, and evaluating these establishes the final result (Chon et al. 2004): 

( 1 /' r \ 2 f 1 a HL+D , o (L+2)! (i-2)! r j < , 

B-WW + 1> ( _' 2 ' 2 o H i " (B6) 

In Section 4.4, we require the evaluation of 
& = $>1)*(2I' + 1)( _\ { L Q )\ (B7) 



(B8) 



This follows along similar lines to Si above. A useful intermediate result is 

4-i(/3) = -4iC3) +csc 2 (/3/2) /^tan^KiO^W, 

Jo 

from which it follows that 

i 5^(21' + l)di'_i(/5) = -<5(cos/3 - 1) + i csc 2 (/3/2). (B9) 
v 

The summation S 2 thus reduces to 

1 / 2 , „ ,„s J , ^ ,L 



S 2 = -y csc 2 (/3/2)di-i(/3)doo(/3)dcos/3, (BIO) 
which, again, vanishes if L ^ Z. Inserting the differential representation of and integrating by parts, we find 
Z I' L V f 1 - ^£±^ for L < Z 



A further result we require in Section 4.4 is 

S 3 = £(-1)*(2Z' + 1)( ^ ^ . (B12) 

Using equations (B2) and (B9), we find 

S 3 = i y C sc 2 (/3/2)d 1 -i(/3)d2-2(/3)dcos/3. (B13) 

In this case, csc 2 (/3/2)g(2_2(/3) i s a polynomial in cos /3 of degree L — 1 and has a root at cos /? = 1. It can thus be expanded as 
a series in with non-zero coefficients for n = 1, . . . ,L — 1, and so the integral in equation (B13) vanishes for L ^ I. The 
summation S3 arises when computing band-power variances from equations (79) and (80); it occurs with prefactors \£lm\ 
and \Blm\ 2 ■ Since we are assuming the weight function is band-limited, L < Z max in S3, so for I 3> Z m ax the contribution from 
S3 vanishes. The final summation we require in Section 4.4 is 

*-E<*+ '>(-> [ -. 0- 

which reduces to 

Si = \ J csc 2 (/V2)4i(/3)d2o(/3)dcos/3. (B15) 
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This time, csc 2 (/3/2)d2o(/3) is a polynomial in cos (5 of degree L — 1 that has a root at cos/3 = —1. The summation Sa can 
thus be expanded in the with n = 1, . . . , L — 1, so that, like S3, it vanishes for L < I. The summation arises from the cross 
term between |3w|iM an d &m in equations (79) and (80), and so this contribution to the band-power variance can also be 
assumed to vanish for I 3> I max ■ 
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